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ABSTRACT 

S' 

Halley's  method  for  the  solution  of  systems  of  equations  is  an  iterative  pro¬ 
cedure  which  converges  cubically  under  favorable  conditions.  The  multivariate  ver¬ 
sion  requires  the  solution  of  two  linear  systems  of  equations  with  the  same  coeffi¬ 
cient  matrix,  following  which  the  correction  vector  is  computed  using  componentwise 
multiplication  and  division  of  vectors.  This  report  describes  a  general-purpose 
computer  program  which  implements  this  method.  The  necessary  first  and  second  de¬ 
rivatives  are  obtained  by  automatic  differentiation,  so  the  user  need  only  supply 
code  defining  the  functions  appearing  in  the  system  of  equations .  The  program  is 
written  in  Pascal-SC,  using  the  new  data  type  HESSIAN  to  represent  dependent  and  in¬ 
dependent  variables.  Numerical  examples  are  given  for  two  simple  systems  of  equa¬ 
tions  to  illustrate  the  use  of  the  program  and  th'.  effectiveness  of  the  method. 
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SIGNIFICANCE  AND  EXPLANATION 


One  of  the  fundamental  problems  of  scientific  computation  is  the  efficient 
numerical  solution  of  systems  of  nonlinear  equations  in  several  variables.  Methods 
are  known  which  converge  rapidly  in  theory,  but  require  first  and  second  partial 
derivatives  of  the  functions  appearing  in  the  system  with  respect  to  the  variables 
involved.  The  necessary  derivatives  can  be  evaluated  automatically,  without  re¬ 
sort  to  numerical  approximations ,  in  programs  produced  by  modem  compilers  which 
permit  user-defined  data  types  and  operators.  Examples  of  such  compilers  are 
Pascal-SC,  Algol  68,  and  ADA  (a  trademark  of  the  U.  S.  Department  of  Defense).  In 
this  study,  Pascal-SC  (Pascal  for  Scientific  Computation)  is  used,  since  it  supports 
accurate  floating-point  arithmetic  for  vectors  and  matrices,  as  well  as  scalars.  The 
method  taken  to  illustrate  these  capabilities  is  Halley's  method,  which  requires 
second  partial  derivatives.  By  use  of  type  HESSIAN,  which  consists  of  the  value  of 
a  function  of  n  variables,  its  gradient  vector  of  first  derivatives,  and  its  Hessian 
matrix  of  second  derivatives  considered  as  a  triple  of  basic  real,  vector,  and  ma¬ 
trix  types,  the  user  need  only  provide  expressions  or  subroutines  for  the  functions 
involved,  and  the  compiler  then  produces  code  for  the  derivatives  needed,  without  re¬ 
sort  to  inaccurate  numerical  or  expensive  symbolic  differentiation.  Since  Halley’s 
method  converges  cubically,  its  speed  can  offsfet  the  overhead  of  calculation  of  the 
necessary  derivatives.  A  general  purpose  programming  system  is  provided  which  con¬ 
sists  of  two  parts:  A  program  which  is  used  to  verify  that  the  coding  of  the  func¬ 
tions  is  correct,  and  another  program  which  solves  the  actual  system.  The  first  pro¬ 
gram  requires  as  input  the  number  of  equations,  the  necessary  operators,  and  expres¬ 
sions  or  subroutines  for  the  functions  involved.  The  second  program  needs  only  the 
number  of  equations  in  the  system  and  the  already  translated  and  verified  code  from 
the  first  program.  A  numerical  example  is  given  to  show  that  this  division  of  labor 
leads  to  efficient  and  effective  numerical  solution  of  a  system  of  nonlinear  equa¬ 
tions  by  the  method  considered. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary 
lies  with  MRC,  and  not  with  the  authors  of  this  report. 


Recession  For 
NTIS  GKA&I  V 

DTIC  TA3 

Unannounced  ^ 


Justification _ 


n 


By- 

Di: 
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1.  Nonlinear  systems  of  equations .  One  of  the  central  problems  of  scien¬ 
tific  computation  is  the  efficient  numerical  solution  of  systems  of  n  equations 


-:t  ion/ 


^-i'.L-iity  Codes 


(1.1)  f .  (x  ,x  , . . . ,x  )  =0,  i  «  1,2,.  ..,n, 

i  1  2  n 

in  n  unknowns  x  ,  x  , . . . ,  x  .  This  is  a  special  case  of  the  opeAatoa  equation 
12  n 

(1.2)  f(x)  »  0, 

in  which  f:D  C  r”  r",  q  €  Rn  denotes  the  zero  vector  0  *  (0,0,..., 0)  ,  and  x  e  Rn 
is  sought,  if  f  is  an  aii-ini  operator, 

(1.3)  f(x)  -  Ax  +  b, 

with  the  matrix  A  *  (a. .)  and  the  vector  b  =  (b, ,b, , . . . ,b  )  given,  then  the  system 
i]  1  2  n 

(1.1)  is  said  to  be  tinea*..  This  important  special  case  is  now  fairly  well  under¬ 
stood  from  a  computational  as  well  as  a  theoretical  standpoint.  Otherwise,  (1.1) 
is  a  nonlinear  system,  and  the  situation  is  quite  different  with  respect  to  theoret¬ 
ical  and  practical  methods  for  solution  than  in  the  linear  case.  Most  of  the  methods 
investigated  to  date  (9),  [10]  involve  some  form  of  iteration,  and  many  also  in¬ 
volve  approximation  of  the  nonlinear  system  by  a  linear  system  during  the  various 
steps  of  the  solution  process .  It  has  been  observed  that  some  solution  procedures 
work  better  than  others  on  a  given  problem,  so  that  in  the  absence  of  a  clear-cut 
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criterion  for  choosing  the  optimal  method,  it  is  advisable  to  have  several  choices 


the  actual  computation j  rather,  the  linear  system 


(2.2)  f'  <xv)av  -  -f  (xV) 

is  solved  for  av,  following  which  the  linear  system 

(2.3)  f  (xu)bu  -  f"(xv)avav 


is  then  solved  for  bu.  Since  the  systems  (2.2)  and  (2.3)  have  the  sane  coefficient 
matrices,  the  decomposition  of  the  Jacobian  matrix  f 1  (xv)  used  to  solve  (2.2)  can 
also  be  used  to  solve  (2.3),  resulting  in  a  saving  of  effort. 

The  total  computation  effort  required  for  one  step  of  Bailey's  method  is  thus: 
1".  Evaluation  of  f(xV),  f ' (xV) ,  f" (xV) ; 

2°.  Solution  of  (2.2)  for  a 

\)  V  \> 

3°.  Evaluation  of  f" (x  )a  a  j 

4°.  Calculation  of  the  Ha&lty  CO VLiCtion  (aV/<aV  +  ^bV)  i 
5°.  Addition  of  the  Halley  correction  to  xv. 

This  sequence  of  operations  is  more  elaborate  than  required  for  Newton's  method, 

which  requires  only  the  evaluation  of  f (xU) ,  f ' (xV) ,  the  solution  of  (2.2)  for  av, 

\)  V  v+1 

and  finally  the  addition  of  a  to  x  to  obtain  x  .  However,  in  favorable  cases, 
the  rate  of  convergence  of  Halley's  method  will  be  cubic,  while  Newton's  method  con¬ 
verges  quadratically.  Thus,  the  greater  effort  required  for  each  step  of  Halley's 
method  could  be  offset  if  fewer  steps  are  required  to  obtain  the  accuracy  desired. 

In  connection  with  the  calculation  of  the  Halley  correction,  an  indeterminate 
form  arises  if  (a'J)i  “  (aV  +  ^t>V)  ^  “  0.  In  this  case,  the  value  of  the  Halley  cor¬ 
rection  is  defined  by  continuity,  and  its  ith  component  is  taken  to  be  0. 

Other  cubically  convergent  iteration  processes  which  use  the  same  information 
as  Halley's  method  are  Chebyihev'i  method  [9],  [10], 


which  is  slightly  less  complicated,  and  the  method  o |j  tangent  hypenbotcu  [9] 


in  which  the  correction  vector  cv  is  obtained  by  solving  the  linear  system 

(2.5)  [f(xv)  +  if"(xv)av]cV  -  -f (xv)  , 
from  which 

(2.6)  xV+1  -  xV  +  cV,  v  «  0,1,2,...  . 

This  method  is  more  complicated  than  Halley's  method  in  that  in  (2.5),  alteration  of 
the  coefficient  matrix  in  (2.2)  is  required.  In  the  scalar  case  (n  »  1) ,  Halley's 
method  and  the  method  of  tangent  hyperbolas  coincide  after  division  of  the  numera¬ 
tor  and  denominator  of  the  Halley  correction  by  aV.  The  abstract  version  of  Halley's 
method  is  defined  in  a  Banach  algebra;  the  Chebyshev  and  tangent  hyperbola  method, 
like  Newton's  method,  do  not  require  multiplication  and  division  of  elements,  and 
hence  can  be  defined  in  the  more  general  setting  of  a  Banach  space  [  2  ) ,  [  9 ) ,  [10) . 

3.  Use  of  automatic  differentiation.  Newton's  method  and  the  methods  of  52 
are  sometimes  shunned  because  it  is  assumed  that  code  has  to  be  supplied  for  the 
derivatives,  or  because  the  functions  f ^  are  defined  by  subroutines,  rather  than 
expressions.  Since  the  rules  for  differentiation  are  well  understood,  however,  the 
computer  can  itself  produce  the  required  code  by  automatic  differentiation  of  the 
given  expressions  or  subroutines  [12].  In  the  case  of  expressions,  programs  capable 
of  obtaining  first  and  second  derivatives  have  been  in  use  for  some  time  [  5 ]  ,  [  7 ]  , 
(lOl •  More  recently,  differentiation  methods  for  subroutines  have  also  been  de¬ 
veloped  [6],  [12],  [14].  Since  the  latter  case  is  the  most  general,  it  will  be 
examined  here. 

In  primitive  computing  languages  such  as  FORTRAN,  automatic  differentiation  re¬ 
quires  interpretation  of  expressions  [5],  [7],  or  precompilation  [6],  [12].  One 
of  the  significant  recent  advances  in  computer  science  has  been  the  development  of 
modern  languages  such  as  Pascal-SC,  in  which  the  performance  of  differentiation  is 
based  on  user-defined  data  types  and  operators  [13],  [14].  To  illustrate  the  basic 
idea,  consider  the  simple  scalar  case  of  a  real  function  f  of  a  single  real  variable 


x.  The  pair  of  values  (f(x),f’(x))  is  the  basic  example  of  a  datum  of  type  GRADIENT 
for  a  given  value  of  x  [13].  Writing  F  -  {f(x),f'(x))  to  represent  an  element  of  this 
new  type  of  data,  the  next  step  is  to  define  the  corresponding  arithmetic  operations 
and  functions  in  a  computable  form.  For  example,  for  G  *  (g(x) ,g* (x) ) ,  addition  and 
multiplication  are  defined  by 

F  +  G  -  (f(x)+g(x),f'(x)+g'(x)), 

(3.1) 

F*G  -  <f(x)*g(x)  ,f(x)*g'  (x)+g(x)*f>  (x) )  , 
respectively.  Similarly,  functions  such  as 

(3.2)  SIN (F)  -  (sin  (f  (x) )  ,f '  (x)  *cos (f (x) ) ) 

are  readily  definable  in  a  form  suitable  for  computational  implementation.  The  in¬ 
dependent  variable  x  is  represented  by  the  GRADIENT  variable  X  ~  (x,l)  ,  and  an  ex¬ 
pression  of  the  form 

(3.3)  F  !-  X*SIN (X  +  4.0)  -  X**3( 

is  then  used  to  obtain  both  the  value  of  the  function  f(x)  «  x  sin  (x  +4.0)  -  x3  and 

2 

its  derivative  f ' (x)  »  x  cos  (x  +  4.0)  +  sin  (x  +  4.0)  -  3x  automatically.  Thus, 
the  user  need  only  supply  the  code  (3.3)  for  the  function  to  be  differentiated,  once 
the  standard  set  of  GRADIENT  operations  and  functions  are  available  [13], 

For  the  present  purpose,  second  derivatives  are  needed,  so  the  type  GRADIENT  is 
extended  to  type  HESSIAN,  a  datum  of  which  is  the  triple  F  -  (f (x) ,f ’ (x) ,f " (x) ) . 

Once  again,  there  is  no  problem  in  the  implementation  of  arithmetic  operations  and 
standard  functions,  for  example, 

F  +  G  -  (f  (x)+g(x)  ,f  <x)+g'  <x)  ,f " (x)+g"  (x) ) 

(3.4) 

F*G  -  <f(x)*g(x)  ,f  (x)  *g'  (x)  +g  (x)  *f '  (x)  ,f  (x)  »g"  (x)+2*f  •  (x)*g'  <x)+g(x)  «f"  (x) )  , 

and 

(3.5)  SIN (F)  «(sin(f  (x) )  ,f'  (x)*cos(f(x)) ,f" (x) *cos (f (x) ) -f • (x)*f’ (x) *sin  (f  (x) ) )  . 

Thus,  given  the  independent  variable  x  as  the  HESSIAN  variable  X  -  (x,l,0),  the 
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evaluation  of  the  expression  (3.3)  yields  the  value  of  the  second  derivative  f"(x) 

m  -  x sin  (x  +  4.0)  +  2  cos  (x  +  4.0)  -  6x  as  well  as  the  values  of  the  function  f(x) 

and  its  first  derivative  f ' (x) . 

Although  the  formulations  of  HESSIAN  operators  and  functions  are  somewhat  com¬ 
plicated,  programming  them  is  no  real  challenge,  and  this  needs  to  be  done  only 
once  and  for  all.  When  available,  these  subroutines  shift  the  burden  of  differen¬ 
tiation  from  the  user  to  the  computing  machine,  which  is  as  it  should  be.  A  com¬ 
plete  package  for  arithmetic  operators  and  standard  functions  has  been  prepared  in 

Pascal-SC  for  the  multivariate  case  in  which  f  is  a  function  of  n  variables,  so 

that  x  -  (x,  ,x_ , . . .  ,x  )  €  Rn,  and  the  HESSIAN  variable  F  is  defined  by 
12  n 

(3.6)  F  =  (f(x) ,7f (x)  ,Hf (x) )  , 
where  Vf(x)  =  f '  (x)  denotes  the  qna.di.lvvt  VictOH 

(3.7)  Vf(x)  =  (3f(x)/3x,  ,3f(x)/3x_,...,3f(x)/3x  ), 

12  n 

and  Hf  (x)  »  f"  (x)  is  the  Huiian  mat/lix 

(3.8)  Hf (x)  =  (32f(x)/3xj3xk) 

of  the  real -valued  function  f  at  x.  The  HESSIAN  variables  X[j)  corresponding  to  the 
independent  variables  x..  are  X[j]  =  (x^  ,e^  ,0)  ,  j  =  l,2,...,n,  where  e^  is  the  jth  , 

unit  vector,  and  0  denotes  the  n*n  zero  matrix.  j 

The  case  of  a  vector-valued  operator  f  is  handled  by  means  of  expressions  or 
subroutines  for  the  n  functions  f ^  appearing  in  (1.1),  defined  to  be  the  correspond¬ 
ing  HESSIAN  variables  F  [i] .  In  this  case,  the  ith  row  of  the  Jacobian  matrix  (1.4) 
is  simply  the  gradient  vector  Vf^(x),  while  the  ith  "panel”  of  the  Hessian  operator 
(1.5)  is  given  by  the  matrix  Hf^x). 

4.  Computation  with  bilinear  operators.  In  the  multivariate  Halley  method 
(2.1),  the  right-hand  side  of  the  linear  system  of  equations  (2.3)  for  bv  is  ob¬ 
tained  by  operating  twice  on  the  vector  av  with  the  bilinear  operator  f"(xv),  where 
the  result  of  the  first  operation  is  a  matrix,  and  the  second  yields  a  vector  [10] . 
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The  way  in  which  HESSIAN  variables  are  defined  makes  it  easy  to  implement  these 
operations.  In  general,  the  bilinear  operator 


(4.1)  B  -  <biJk> 

will  be  considered  to  be  composed  of  n  matrices 

<4‘2>  B1  “  (bl jk> '  B2  ■  <b2jk> .  Bn  ‘  (bnjk)' 

which  will  be  called  i-paneti,  or  simply  paneti  of  B.  For  a  vector  x£  r”,  the 
matrix 

(4.3)  A  -  (aAj)  -  BX  -  ) 

T 

will  have  rows  A^  given  by  the  matrix-vector  product 

(4.4)  A^  -  B^x,  i  »  1,2,... ,n 

Once  the  matrix  A  is  formed  by  computing  the  vectors  (4.4)  ,  then  the  vector 

n  n 

(4.5)  y  -  Ax  -  Bxx  -  (  £  l  b.-X-x,  ) 

j-1  k-1  13*  *  3 

is  obtained  by  a  single  additional  matrix-vector  multiplication.  Here,  Bi  •  Hfi(xu), 

(4.6)  -  Hfi(xV)aV,  i  -  1,2,. ...n, 
and  thus 

(4.7)  f  (x  )a  a  °  Aa  , 

so  the  required  vector  is  obtained  by  a  total  of  n  +  1  matrix-vector  multiplications 
In  Pascal-SC,  vectors  and  matrices  are  stored  row-wise,  so  no  transposition  is 
required  when  forming  the  matrix  A  from  the  vectors  A^  given  by  (4.6)  [IS] • 

5.  Programming  Halley's  method  in  Pascal-SC.  Central  to  Pascal-SC,  as  well  as 
to  Pascal  [  1  J ,  is  the  concept  of  a  data.  type..  In  Pascal-SC,  vectors  and  matrices 
over  type  REAL  (the  set  of  floating-point  numbers)  are  considered  to  be  the  standard 
types  RVECTOR  and  RMATRIX,  respectively  (15).  Following  the  conventions  of  Pascal-SC 
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P 


n -dimensional  vectors  and  matrices  are  declared  in  the  program  heading  by: 


CONST  DIM  =  n; 

TYPE  DIMTYPE  =  I ..DIM; 

(5-1)  RVECTOR  =  ARRAY[DIMTYPE]OF  REAL; 

R MATRIX  =  ARRAY[0!MTYPE]0F  RVECTOR; 

and  a  number  of  ordinary  operations  of  matrix  and  vector  algebra  are  implemented 
[8  1.  Following  (5.1),  type  HESSIAN  is  declared  by: 

(5.2)  TYPE  HESSIAN  =  RECORD  F:  REAL  ;DF :  RVECTOR ;HF:  R MATRIX  END; 

so  that  if  an  expression  or  the  result  of  a  subroutine  for  computing  f(x)  is  as¬ 
signed  to  the  hessian  variable  F,  one  has 

(5.3)  F.F  =  f(x),  F.DF  =  Vf(x),  F.HF  =  Hf(x). 


Step  1°  of  Halley's  method  as  outlined  in  §2  is  thus  taken  care  of  simply  by 

expressing  the  independent  variables  x,  ,  x,,...,  x  and  the  values  of  the  functions 

l  ^  n 

fl'  f 2  ' '  "  '  fn  in  the  sYsteo  IT-1)  as  HESSIAN  variables.  For  example,  consider  the 
simple  system  of  equations 

e'x+y  -  o.l  *  0, 

(5.4) 

e~x~y  _  0<1  _  Qf 

investigated  by  Cuyt  and  Van  der  Cruyssen  [ 2  ] ,  [ 4  I .  The  variables  involved  would 
be  declared  to  be  HESSIAN  in  the  heading  of  the  program  (see  Appendix  A)  by 

(5.5)  VAR  X.Y.F.G:  HESSIAN; 


and  the  functions  corresponding  to  the  left-hand  sides  of  (5.4)  by 


(5.6) 


F  :*  HEXP(-X+Y)  -  0.1 ; 
G  ;=  HEXP(-X-Y)  -  0.1; 


in  the  body  of  the  program.  (The  user  is  free  to  name  and  order  both  the  independent 
and  dependent  variables  in  any  convenient  manner.  A  more  systematic  approach  will  be 


-  8  - 


discussed  in  the  next  section.)  In  (5.6),  the  convention  that  the  names  of  standard 
functions  for  type  HESSIAN  begin  with  “H“  has  been  followed.  Evaluations  of  the 
expressions  in  (5.6)  requires  the  following  HESSIAN  operators  and  functions: 

OPERATOR  -  (H:  HESSIAN)  RES:  HESSIAN; 

OPERATOR  -  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

(5.7)  OPERATOR  -  (H:  HESSIAN;R:  REAL)  RES:  HESSIAN; 

OPERATOR  +  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

FUNCTION  HEXP(H:  HESSIAN):  HESSIAN; 

source  code  for  these  is  given  in  the  heading  of  the  program  listed  in  Appendix  A. 
The  first  operator  calculates  -H,  the  second  HA-HB,  and  so  on. 

The  independent  variables  X,Y  are  initialized  as  follows:  Their  function  value 
(or  simply  value)  parts  x.F  and  Y.F  are  given  initially  by  input  from  the  user,  and 
are  subsequently  calculated  by  the  Halley  iteration.  Their  gfia.die.nt  parts  x.DF  and 
Y.DF  are  assigned  the  constant  values 

(5.8)  X.DF[1  ]:»1 ;  X.DF[2]:=0;  Y.0F[1 ]:=0;  Y.DF[2]:=1; 
as  are  their  He&itan  parts 

(5.9)  X.HF:=MRNULL ;  Y.HF:*MRNULL; 

where  MRNULL  is  a  standard  Pascal-SC  function  which  returns  the  zero  matrix. 

Execution  of  the  statements  (5.6)  thus  completes  step  1°  of  Halley's  method. 

For  step  2°,  the  Jacobian  matrix  of  the  system  (5.4)  is  needed.  It  is  assumed  that 
the  declaration 

(5.10)  VAR  JAC.L.M:  RMATRIX ;  A.B ,V:  RVECTOR;  NRS :  BOOLEAN; 

is  in  the  heading  of  the  program,  where  JAC  denotes  the  desired  Jacobian.  It  is 
obtained  by  the  assignments 

(5.11)  JAC[1]:=F.DF;  JAC[2]:=G.DF; 

since  its  rows  are  the  gradient  vectors  Vf(x)  and  7g(x)  ,  respectively. 


Similarly,  letting  V  denote  the  right-hand  side  of  the  linear  system  (2.2),  one 


has 

(5.12)  V[l]:-  -F.F;  V[2]:=  -G.F ; 

v 

and  all  that  remains  is  to  solve  (2.2)  for  A  =  a  by  means  of  a  standard  procedure 
for  solving  linear  systems  of  equations,  such  as 

(5.13)  SOLVLN(DIM,JAC,M,V,A,NRS); 

in  which  the  decomposition  of  JAC  is  stored  as  the  matrix  M,  and  additional  right- 

hand  sides  will  be  expected  as  long  as  NRS  =  FALSE.  This  completes  step  2°. 

v  v 

In  step  3°,  the  matrix  f “ (x  )a  is  computed  as  the  matrix  L.  This  is  accom¬ 
plished  by  means  of  the  assignments 

(5.14)  L[1 ] :=  F.HF*A;  l[2]:=  G.HF*A; 

using  the  standard  Pascal-sC  operator  *  for  matrix  by  vector  multiplication  115). 

V  V  V 

With  this  result,  the  vector  V  =  f"(x  )a  a  is  obtained  from 

(5.15)  V:=  L*A; 
thus  completing  step  3°. 

The  calculation  of  the  Halley  correction  (step  4°)  requires  first  the  compu- 
tation  of  the  vector  B  =  b  by 

(5.16)  SOLVLN(DIM, JAC ,M,V,B,NRS); 

where  V  is  now  obtained  from  (5.15),  and  NRS  =  TRUE.  Addition  of  vectors  and 
multiplication  of  vectors  by  real  numbers  are  standard  in  Pascal-SC;  however , 
componentwise  multiplication  and  divison  of  vectors  are  not.  Thus,  the  corres- 
jjonding  operators  *,/  must  be  defined  in  the  heading  of  the  program  by 

OPERATOR  *  (VA.VB:  RVECTOR)  RES:  RVECTOR; 

VAR  U:  RVECTOR; I :  DIMTYPE; 

<5-17)  BEGIN  FOR  I :=1  TO  OIM  DO  U[I]:=VA[I]*VB[I]; 

RES : *  U 

END; 


ft 
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and 


OPERATOR  t  (VA.VB:  RVECTOR)  RES:  RVECTOR ; 

VAR  U:  RVECTOR;  I:  DIMTYPE ; 

BEGIN  FOR  I:-1  TO  DIM  DO 

(5.18)  IF  (VA[I]  =  0)  AND  (VBtl]  =  0)  THEN  U[I]:*0 

ELSE  U[I]:»  VA[1]/VB[I]; 

RES:*  U 

ENO; 

respectively.  The  division  operator  is  tailored  to  yield  0  as  the  limit  of  the 

v 

Halley  correction  as  (a  )  -*  0,  which  is  valid  as  long  as  there  is  a  neighborhood 

v 

of  a  which  does  not  contain  points  for  which  the  denominator  of  the  Halley  correc¬ 
tion  is  zero  while  the  numerator  is  nonzero.  By  the  use  of  the  operators  (5.17) 
and  (5.18),  the  Halley  correction  is  the  vector  V  given  by 

(5.19)  V:=  (A*A)/(A  +  0.5*8); 

which  completes  step  4°. 

The  final  step  of  one  Halley  iteration  is  then 

(5.20)  X . F : *  X.F  +  V[1 ] ;  Y.F:«  V.F  +  V[2j; 

after  which  another  iteration  can  be  performed,  if  desired.  A  set  of  typical  numer¬ 
ical  results  for  this  problem,  using  the  program  in  Appendix  A,  is  given  in  Appendix 
B. 

6.  A  more  general  approach;  Type  SYSTEM.  in  the  simple  example  discussed 
in  the  previous  section,  it  was  convenient  to  use  the  ordinary  notation  X,Y  for 
the  independent  variables  involved,  and  F,G  for  the  dependent  variables  corres¬ 
ponding  to  the  system  (5.4).  For  larger  systems,  it  is  helpful  to  adopt  a  more 
formal  notational  convention,  in  terms  of  which  a  general-purpose  program  can  be 
developed.  This  is  done  by  the  introduction  of  the  data  type  SYSTEM,  which  is 
declared  by: 

(6.1)  TYPE  SYSTEM  *  ARRAY[DIMTYPE]OF  HESSIAN; 
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by  the  use  of  this  data  type,  a  declaration  of  the  form 

(6.2)  VAR  X.F:  SYSTEM; 

can  be  used,  for  example,  to  introduce  independent  variables  X (1) ,X (2) , . . . ,x (DIM] , 
and  dependent  variables  F [1] ,F [2] , . . . ,F [DIM] .  Thus,  the  system  of  equations 

4  4  4 

16x^  +  i6x^  +  x^  -  16  =  0, 

2  2  2 

(6.3)  Xj^  +  x2  +  x3  -  3  =  0, 


taken  from  [13]  would  be  coded  as 

F[1 ]:-  1 6*( X[1 ]**4 )  +  16*(X[2]**4)  +  X[3]**4  -  16; 

(6.4)  F[2]:=  X[1 ]**2  +  X[2]**2  +  X[3]**2  -  3; 

F[3]:-  X [ 1 ]**3  -  X[2 ] ; 

(parentheses  are  necessary  in  (6.4),  since  •*  and  *  have  the  same  priority  in  Pascal-SC 
(11).  Once  the  statements  (6.4)  are  executed,  the  value  F[i).F  of  each  function, 
the  ith  row  of  the  Jacobian  matrix  of  the  system  F  [  i ] . DF ,  and  the  ith  panel  F[i).HF 
of  the  Hessian  operator  are  all  at  the  disposal  of  the  programmer. 

7.  The  programs  SYSTEST  and  HALSYS.  A  brief  description  of  two  Pascal-SC 
programs  w.iich  can  be  used  to  investigate  the  .ipplicatron  of  Halley's  method  to 
systems  of  equations  will  now  be  presented.  First  of  all,  good  programming  prac¬ 
tice  requires  that  the  expressions  and  subroutines  provided  by  the  user  produce  the 
correct  values  of  the  functions  f^(x)  appearing  in  (1.1).  The  program  SYSTEST  is 
provided  for  this  purpose.  To  use  this  program,  the  file  SYSTEST. S  containing  its 
source  code  (see  Appendix  C)  is  created  from  SYSTEST. PROG  and  edited  to  include  the 
correct  value  of  DIM  and  expressions  or  subroutines  for  the  functions  F  [I ] .  Source 
code  for  the  required  HESSIAN  operators  and  functions  is  obtained  from  the  file 
HESS_PAKET,  which  contains  the  22  arithmetic  operators,  5  power  operators,  and  the 
functions  HABS,  HSQRT,  HEXP,  HU),  HARCTAN,  HSIN,  and  HCOS  [14).  As  explained  in 
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[14] >  the  user  may  add  any  other  needed  HESSIAN  functions,  procedures,  or  operators 
to  the  program  heading  if  needed  to  supplement  the  standard  ones.  For  example,  the 
system  (6.4)  requires  that  source  code  for  the  following  operators  appears  in  the 
program  heading: 

OPERATOR  *  (K:  INTEGER;H:  HESSIAN)  RES:  HESSIAN; 

OPERATOR  **  (R:  REAL ;K:  INTEGER)  RES:  REAL; 

OPERATOR  **  (H:  HESSIAN;K:  INTEGER)  RES:  HESSIAN; 

(7.1) 

OPERATOR  +  {HA.HB:  HESSIAN)  RES:  HESSIAN; 

OPERATOR  -  (H:  HESSIAN-.K:  INTEGER)  RES:  HESSIAN; 

OPERATOR  -  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

and  one  sets  DIM  ■  3.  After  this  program  is  compiled,  upon  execution  it  will  re¬ 
quest  initial  values  for  X(l]  .F, . . .  ,X  [DIM]  .F,  and  then  print  out  the  values  F[I).F 
of  the  functions  in  the  system,  the  values  of  the  first  derivatives  DF[I]/DX[J)  = 
F[I].DF[J],  and  the  second  derivatives  D2P[I]/DX[J]DX[K]  «  F[IJ  ,HF[J,K] .  A  typical 
set  of  output  for  the  system  (6.4)  is  given  in  Appendix  D. 

Once  the  correctness  of  coding  for  the  system  has  been  verified,  the  translated 
code  for  the  procedure  FCOMP(VAR  X,F:  SYSTEM; DIM:  INTEGER)  is  entered  in  the  external 
library  HALLE  Y_LIB  as  subroutine  number  777  [15]  from  the  intermediate  code  (20  for 
the  program  SYSTEST.  HALLEY_LIB  also  contains  pretranslated  code  for  the  linear 
equation  solver  SOLVTA  (subroutine  number  776)  and  the  standard  functions  of  Pascal-SC 
(15J.  The  source  code  file  HALSYS.S  is  created  from  HALSYS . PROG  with  DIM  set  to  its 
correct  value,  and  the  program  HALSYS  compiled  with  reference  to  HALLEYJLIB  to  bring 
in  the  code  for  the  system  being  solved.  The  text  of  HALSYS. PROG  is  given  in  Appen¬ 
dix  E. 

The  program  HALSYS  carries  out  the  actual  Halley  iteration.  First,  it  asks  for 
initial  values  of  X[l) ,F, . . . ,X [DIM] .F,  and  then  prints  these  and  the  corresponding 
function  values  F [1] .F, . . . ,F [DIM] .F.  The  user  then  receives  the  query: 
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(7.2) 


RESTART  (R)  OR  ITERATE  (Y/N)? 


The  response  "R"  will  result  in  the  request  for  another  set  of  initial  values,  "Y" 
will  give  the  results  of  one  Halley  iteration,  while  "N"  will  terminate  the  program 
and  return  control  to  the  operating  system.  After  each  step  of  Halley's  method, 
the  query  (7.2)  will  be  sent  to  the  user.  Of  course,  the  user  can  introduce  a  more 
automatic  method  for  controlling  the  iteration  by  editing  the  file  HALSYS.S  before 
compilation.  Typical  results  for  the  system  (6.3)  are  given  in  Appendix  F.  It 
is  interesting  to  compare  these  with  the  corresponding  results  obtained  by  Newton's 
method,  and  given  in  (13).  In  the  latter  case,  eight  iterations  were  required  to 
reduce  the  function  values  to  zero,  as  compared  to  6  by  Halley's  method.  However, 
a  few  iterations  are  spent  in  each  calculation  chasing  a  small  roundoff  error  in 
the  function  values;  the  function  values  are  actually  negligible  after  5  iterations 
of  Newton's  method  and  3  iterations  of  Halley's  method. 

6.  Implementation  details.  The  programs  described  in  this  report  were  created 
and  tested  using  the  Paacal-SC  compiler  developed  at  the  University  of  Karlsruhe  for 
the  Zilog  MCZ-1  microcomputer  using  the  RIO  2.06  operating  system  of  zilog,  Inc.  No 
other  claims  of  correctness  or  usability  are  made. 
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APPENDIX  A 


A  Pascal-SC  program  for  the  solution  of  system  (S.4)  by  Halley's  method. 


PROGRAM  HALLEY  (INPUT,  OUTPUT)* 

CONST  DIM  -  2* 

TYPE  DIMTYPE  -  1..DIM) 

R VECTOR  -  ARRAY  [DIMTYPE]  OP  REAL* 

RMATRIX  -  ARRAY  [ DIMTYPE]  OP  RVECTOR* 

HESSIAN  “  RECORD  Pi  REAL*  DP  I  RVECTORjHFi  RMATRIX  END* 

VAR  X,Y,P,Gi  HESSIAN |Ci  CHAR*V,A,B:  RVBCTORl JAC,L,M:  RMATRIX* 

NRSi  BOOLEAN* 

(*  The  following  are  standard  Pascal-SC  matrix  and  vector  functions  and 
operators  from  MR_PAKET.  •) 

FUNCTION  MRNULLi  RMATRIX*  (*  This  returns  the  aero  matrix.  *) 

VAR  X,Jt  DIMTYPE* 

Ci  RMATRIX* 

BEGIN 

FOR  Ii-1  TO  DIM  00 
FOR  Ji-1  TO  DIM  DO 
C(I,J]  J-  0* 

MRNULL  i-  C 
END* 

OPERATOR  +  (A,B|  RVECTOR)  RESi  RVECTOR* 

VAR  It  DIMTYPE* 

BEGIN  FOR  Ii-1  TO  DIM  DO  A[I]  i-  A[I]+B[IJ» 

RES  i-  A 
END* 

OPERATOR  •  (At  REAL*  Bi  RVECTOR)  RESi  RVECTOR* 

VAR  It  DIMTYPE* 

BEGIN  FOR  It-1  TO  DIM  DO  B(I]  !-  A*B[I]| 

RES  i-  B 
END* 

OPERATOR  *  (At  RMATRIX*  Bi  RVECTOR)  RESi  RVECTOR* 

VAR  It  DIMTYPE* 

BVARi  RVECTOR* 

BEGIN 
BVAR  l-  B* 

FOR  Ii-1  TO  DIM  DO 

B[I]  t-  SCALP  ( A [ I ] , BVAR, 0 ) * 

RES  t-  B 
END* 
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(*  Special  operators  for  componentwise  multiplication  and  division  of  vectors, 
see  (5.17)  and  (5.18).  •) 

OPERATOR  *  (VA.VBi  RVECTOR)  RES  I  RVECTORj 
VAR  U:  RVECTOR;  1 1  DIMTYPE) 

BEGIN  FOR  Is»1  TO  DIM  DO  Utl]  I-VAtU  *VB[I)  ; 

RES i “U 

END; 

OPERATOR  /  (VA,VBf  RVECTOR)  RES:  RVECTOR; 

VAR  Ut  RVECTOR;  It  DIMTYPE; 

BEGIN  FOR  I:-1  TO  DIM  DO 

IF  (VA[I]-0)  AND  (VBtU-O)  THEN  U[I):-0 
ELSE  D(U:-VA(I]/VB[I1; 

RES.'-O 

END; 

(*  The  required  HESSIAN  operators  and  function  (5.7)  for  the  evaluation  of 
the  system  (5.6)  follow.  *) 

OPERATOR  +  (HA,HB:  HESSIAN)  RESt  HESSIAN; 

VAR  I ,  J:  DIM.*PE;Ui  HESSIAN; 

BEGIN  U.Ft-HA.F+HB.FlFOR  It-1  TO  DIM  DO 
BEGIN  U.DF(I) t “HA.DF [I] +HB.DF [I] ; 

FOR  Jt«»  TO  DIM  00 

O.HFII]  (J)t-HA.HF(I)  [J]+HB.HF[I]  [J] 

END; 

RES:-U 

END; 

OPERATOR  -  (Hi  HESSIAN)  RES:  HESSIAN; 

VAR  I,J:  DIMTYPE ; 0 :  HESSIAN; 

BEGIN  U.F:—  H.FiFOR  I:»1  TO  DIM  DO 
BEGIN  U.DF(I):— H.DFtI]  ; 

FOR  Jj-1  TO  DIM  DO 
U.HFU]  (J]t— H.HFlI]  (J] 

END; 

RES:"U 

END; 

OPERATOR  -  (Ht  HESSIAN  ;R:  REAL)  RES  I  HESSIAN; 

VAR  Ut  HESSIAN; 

BEGIN  U.F:-H.F-R;U.DF:-H.DF;U.HFi-H.HF; 

RESt-0 


OPERATOR  -  (HA.HBi  HESSIAN)  RESi  HESSIAN ; 

VAR  I,Jt  DIMTYPE)Ut  HESSIAN) 

BEGIN  U  •  F  x “HA •  F-HB •  P J FOR  Ii-I  TO  DIM  DO 
BEGIN  U.DF(I];«HA.DF[I]-HB.DF[Il; 

FOR  J|-1  TO  DIM  DO 

U.HF[I]  U]l-HA.HF[I]  [J]-HB.HFtU  [J] 


END) 


RES:-U 


END) 


FUNCTION  HEXP { H s  HESSIAN))  HESSIAN; 

VAR  I,J)  DIMTYPE)U>  HESSIAN) 

BEGIN  U.Fi~EXP(H.F); 

FOR  Il-1  TO  DIM  DO 

BEGIN  U.DFlI] x»U.F*H.DFlI)) 

FOR  Jl-1  TO  I  DO 

BEGIN  U.HF[I]  [J]  x-U.F*H.HF[I]  [J]+U.DF[I]  *H.DF(J]  ) 

IF  IOJ  THEN  U.HF[J]  (I)  i-U.HF[IJ  [Jl 

END; 

END; 

HEXP : -U 
END; 

(*  The  next  procedure  eolvea  a  linear  system  of  equations  with  coefficient 
matrix  JAC  and  right-hand  side  V  for  the  solution  vector  S.  The  decaeposi- 
tion  of  JAC  is  stored  as  the  matrix  M.  If  NRS  “  FALSE,  then  additional 
right-hand  aides  are  expected.  Fre-translated  code  for  this  procedure  is 
stored  in  the  library  HALLEY_LIB  as  subroutine  number  776.  * ) 

PROCEDURE  SOLVLN  (DIM s INTEGER; VAR  JAC,MlRMATRIX;VAR  V,St  RVBCTORlNRS x BOOLEAN )  ; 
EXTERNAL  776) 

BEGIN  (*  Initialisation  of  gradients  and  Hessians  of  independent  variables.  *) 

X.DP[1]:-1)X.DF[2]i-0)Y.DF(1]:-0)Y.DF[2] t-1; 

X.HFi“MRNULL)Y.HF:^!RNULL) 

C:*'R® ) WHILE  C  -  'R'  DO 

BEGIN  (*  VALUE  INITIALIZATION  *) 

WRITELN ( 'ENTER  X, Y 1 ) I  READ (X. F, Y. F ) ) 

Cl-'Q';  WHILE  C  -  '2'  DO 

BEGIN  (•  MAIN  PROGRAM  •) 

(*  Calculate  function  values  and  derivatives.  Print  values  of 
independent  and  dependent  variables.  *) 

Fx-HEXP(-X+Y)-0. 1; 

G:»HEXP(-X-Y)-0.1; 
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WRITELN('(X,Y)  -  ('fX.F,*  ,  ,,Y.F,*),)» 

WRXTELN  ( '  (F#G)  -  <’,F.F,'  ,  *,G.F,,),)| 

WRITELNt  'RESTART  (R)  OR  ITERATE  (Y/N  )  ?  *  )  |  READ(C,C ) > 

(*  In  response  to  "H",  the  program  will  ask  for  new  initial  values  of 
the  independent  variables)  an  input  of  *Y*  will  result  in  one  Halley 
iteration  being  performed,  while  *N"  will  terminate  the  program  and 
return  control  to  the  operating  system.  Two  characters  are  read, 
since  the  first  read  from  the  console  will  always  be  a  blank,  the 
second  being  the  character  entered  by  the  user  in  response  to  the 
prompt  *) 


WHILE  C  -  'Y '  DO 

BEGIN  (*  HALLEY  ITERATION  *) 

(*  Construct  Jacobian  matrix  JAC  and  right-hand  side  V  of  (2.2).  * ) 

JAC[1]:-F.DFjJAC[2]:-G.DF>V[1) )—  F.F|V[2Js— G.Fi 

(•  Solve  (2.2)  for  A.  *) 

NRS  >  -FALSE I SOLVLN  (DIM,  JAC  ,M,  V,  A  ,NRS  )  ) 

(*  Construct  the  right-hand  side  Vof  (2.3).  *) 

L (1 )  i -F.HF*A|L  [2]  s^3.HF*A)V:-L*A) 

(*  Solve  (2.3)  for  B.  •) 

NRS)— TRUE) SOLVLN (DIM, JAC, M,V,B, NRS) > 

(*  Compute  the  Halley  correction  V  and  update  independent  variables.  *) 

Vi-A*A/(A+0.5»B)l 
X.F)-X.F+V (1 J ) Y.F:— Y.F+V [2] ) 

Ct-'Q* 

END)  (•  HALLEY  ITERATION  •) 

END)  <»  MAIN  PROGRAM  *) 


END  (*  VALUE  INITIALIZATION  *) 


APPENDIX  B 

Output  of  the  program  in  Appendix  A  tor  tha  initial 
valuaa  X.p  -  4.3,  Y.r  -  2.0. 

INITIAL  VALUES 

(X,Y)  -  (  4. 30000000000E+00  ,  2.  OOOOOOOOOOOR+OO ) 

(F,G)  -  {  2.588437230001-04  ,  -9.816369522301-02) 

RESULTS  OF  ITERATION  NUMBER  1 

(X,T>  -  (  3. 336 1 S528246E+00  ,  1.03597241993E+00 ) 

(F,G)  -  (  2.4051 1813000K-04  ,  -8.73756488903E-02) 

RESULTS  OF  ITERATION  NUMBER  2 

(E.Y)  -  (  2. 5608180 0937E+00  ,  2. S9679794972E-01 ) 

(F,G)  -  (  1.44792583000E-04  ,  -4.042372200388-02) 

RESULTS  OF  ITERATION  NUMBER  3 

(E.Y)  -  (  2.3081 7563469E+00  ,  5.68378530700E-03 ) 

(F,G)  -  (  9. 32479600000B-06  ,  -1.121100995708-03) 

RESULTS  OF  ITERATION  NUMBER  4 

(X,Y)  -  (  2.302585151 19E+00  ,  6. 12025800000E-08) 

(F,G)  -  (  3. OOOOOOOOOOOE-IO  ,  -1. 19396000000E-08) 

RESULTS  OF  ITERATION  NUMBER  5 

(X,Y )  -  (  2 . 30258509299S+00  ,  4.57643B40C00B-12) 

(F,G)  -  (  O.OOOOOOOOOOOE+OO  ,  0.000000000008+00) 


APPENDIX  C 


The  program  SYSTEST  for  testing  correctness  of 
HESSIAN  systems  of  equations. 


The  source  code  for  this  program  is  in  the  file  SYSTEST. PROG. 


PROGRAM  SYSTEST(INPUT, OUTPUT) ) 

CONST  DIM  =*  #;  (*  Replace  "#“  by  the  dimension  of  the  system  teBted.  *) 

TYPE  D1MTYPE  -  1..DIM) 

RVECTOR  -  ARRAY [DIMTYPEjOF  REAL) 

RMATRIX  -  ARRAY tDIMTYPE) OF  RVECTOR) 

HESSIAN  -  RECORD  F:REAL)DF:  RVECTOR /HF:  RMATRIX  END; 

SYSTEM  -  ARRAY [DIMTYPEJOF  HESSIAN) 

VAR  X.Fi  SYSTEM)  JAC:  RMATRIX)I,J,K:  DIMTYPE)C:  CHAR) 

PROCEDURE  FCOMP (VAR  X,F:  SYSTEM) DIM:  DIMTYPE ) ) 

l*  Insert  source  code  for  the  operators,  functions,  and  procedures  required 
for  computation  of  the  system  being  tested  here,  for  example,  the  source  code 
for  the  operators  (7.1)  in  the  case  of  the  system  (6.4).  Source  code  for 
HESSIAN  operators  and  functions  is  in  the  file  HBSS_PAKET.  *) 

BEGIN  (*  SYSTEM  DEFINITION  •) 

(•  Insert  code  defining  the  system  to  be  tested  here,  for  example,  the  system 
(6.4): 

F  [1) :“16*(X [1 )**4)+16*(X [2) **4)+X [3] **4-16) 

F  [2  J  :-X  [11  ««2+X  [2]  **2+X  [3]  **2-3) 

F[3J  :-X[11**3-X[2J  )  •) 

END)  (*  SYSTEM  DEFINITION  *) 

(*  The  following  are  standard  Pascal*SC  functions  from  MR_PAKET:  MRNUI.L 
returns  the  zero  matrix,  and  MRID  returns  the  identity  matrix.  These  are 
used  for  initialization  of  the  Hessian  and  gradient  parts  of  the 
independent  variables,  respectively.  *) 
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FUNCTION  MRNULLx  RMATRIX; 

var  i,j«  dimtype; 

Cx  RMATRIX ; 

BEGIN 

FOR  Ij-1  TO  DIM  DO 
FOR  J:»1  TO  DIM  DO 
C(I,J]  I-  0; 

MRNULL  !-  C 
END; 

FUNCTION  MRXD:  RMATRIX; 

VAR  1,J s  DIMTYPE; 

C;  RMATRIX; 

BEGIN 

FOR  Ix-1  TO  DIM  DO 
FOR  J:-1  TO  DIM  DO 
IF  I-J  THEN  C[I,J1  X-  1 
ELSE  C(I,J]  X-  0; 

MRID  X-  C 
END; 


begin  (*  Initialization  of  gradient*  and  Haaslana  of  independent  variable*.  *) 
JAC:-MRID;FOR  Ix-1  TO  DIM  DO 
BEGIN  X[I].OF,-JAC[I]»X(I).HFx-HRNULL 
END; 

Cx— 'Y* ; WHILE  C  -  'Y'  DO 
BEGIN  (*  MAIN  PROGRAM  *) 


(*  Input  portion  *) 

writeln Center  independent  variables'); 

FOR  Ix-1  TO  DIM  DO  READlX  {I )  . P); 

(*  Calculation  of  function  values  and  derivatives.  *) 
FCOMP(X,F,DIM);WRITELN( '  ')» 

(*  Output  Portion  •) 

FOR  Ix-1  TO  DIM  DO 

BEGIN  WRITELN ('FUNCTION  VALUE: ’ );  WRITELN (*Ft,,I»2,,l  •  ',FlIJ.F)l 
WRITELN ('FIRST  DERIVATIVES: ' );FOR  Jx-  1  TO  DIM  DO 
WRITELN ( 'DF [',1x2,'] /DX [ ' , Jx  2, ' ]  -  ■ ,F (I J .DF tJl  ) ; 

WRITELN ( 'SECOND  DERIVATIVES: ' );FOR  J:-1  TO  DIM  DO 
FOR  Kx-1  TO  J  DO 

WRITELNCDZFlMxa.'l/DXl'rJxI/lDXl'.KxJ,']  -  '  ,F  [I)  .HF  [J,K]  ) 
END;  (*  Output  Portion  *) 


23  * 


WRITELN  (  ' DO  YOU  WANT  TO  ENTER  MORE  VALUES  ( Y/N  )  ?  •  )  ;READ(C,C) 


r 


END;  (*  MAIN  PROGRAM  *) 


END. 


EXAMPLE:  Source  code  for  the  HESSIAN  operators  (7.1)  for  the  system  (6.4). 


OPERATOR  *  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

VAR  I,J:  DIMTYPE/U:  HESSIAN; 

BEGIN  U . F : “HA • P+HB • F ; FOR  I:-1  TO  DIM  DO 
BEGIN  U.DFtI]  :-HA.DF(I]  «-HB.DF[I]  ; 

FOR  J:-1  TO  DIM  DO 

U.HFtl.J] :-HA.HFtI,J]+HB.HF[I,J] 

END; 

RES : “U 

END; 

OPERATOR  -  (H:  HESSIAN ;K:  INTEGER)  RES:  HESSIAN; 
VAR  U:  HESSIAN; 

BEGIN  U.F:-H.F-K;U.DF:-H.DF;U.HF:-H.HF; 

RES:-U 

END; 

OPERATOR  -  (HA.HB:  HESSIAN)  RES:  HESSIAN; 

VAR  I.J:  DIMTYPE;U:  HESSIAN; 

BEGIN  U.F:-HA.F-HB.F;FOR  t:-1  TO  DIM  DO 
BEGIN  U.DF[I] :-HA.DF(I]-HB.DF(I] ; 

FOR  Jim1  TO  DIM  DO 

U.HF[I,J] :-HA.HF[I,J]-HB.HF[I,J) 

END; 

RES:-U 

END; 

OPERATOR  *  (K:  INTEGER; H :  HESSIAN)  RES:  HESSIAN; 
VAR  I.J:  DIMTYPE;U:  HESSIAN; 

BEGIN  U.F:-K*H.F;FOR  Ii-I  TO  DIM  DO 
BEGIN  U.DF[I] :-K*H.DFtI] ; 

FOR  J:“1  TO  DIM  DO 
U.HFtl.J] :»K*H.HFtI,J] 

END; 

RES : -U 
END; 
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OPERATOR  ••  (Ri  REALlXl  INTEGER)  RESl  REAL  I 
VAR  Li  INTEGER; Ut  REAL; 

BEGIN  IMt  <■  0  THEN  Ui-I/R; 

IE  (K  -  0)  OR  (R  -  1)  THEN  Ui-1 
ELSE  IF  K  -  1  THEN  Ui-R 

ELSE  BEGIN  Li  -ABS (X )  ;Ut  “1  (REPEAT  IF  L  MOD  2  -  1 
THEN  Ui-R*U;Li-  L  DIV  2;  IF  L  <>  0 
THEN  Ri-R*R  UNTIL  L  -  0; 

IF  X  <  0  THEN  Ut-1/U 

END; 

RESi-U 

END; 

OPERATOR  •*  (Hi  HESSIAN iKl  INTEGER)  RESl  HESSIAN; 

VAR  I,Ji  DIMTYPS;H,MMiREAL;Ut  HESSIAN; 

BEGIN 

IF  K  -  0  THEN 
BEGIN  IF  H.F  •  0  THEN 

BEGIN  WRITEWI 'EXPONENTIATION  ERROR  0**0 '  )  ;SVR(0  ) 

END 

ELSE  BEGIN  U.Fl-1;FOR  Il-I  TO  DIM  DO 

BEGIN  U.DF(Ili-0;FOR  Jl-1  TO  DIM  DO  U.HF(I,J)t-0 
END; 

END; 

END 

ELSE  IF  K  -  1  THEN  Ul-H 
ELSE  IF  X  -  2  THEN 

BEGIN  U.Fi-H.F*H.F;Mi-2*H.F|FOR  Ii-I  TO  DIM  DO 
BEGIN  U.OF[I]|-M*H.OF[Z];FOR  Jl-1  TO  I  DO 

BEGIN  U.HF[I,J]|-M*H.HF[I,J)+2*H.DF[I]*H.DF[J]| 
IF  I  <>  J  THEN  U.HF[J,I]s-U.HF[I,J] 

END; 

END; 

END 

ELSE  BEGIN  MNi-H.F**<K-2 ) lMi-H.F»MM;U.Ft-H.F*M; 

Ml-K»M|MMt-X*<K-1)*MNlFOR  Il-I  TO  DIM  DO 
BEGIN  U.DF[I]i-M*H,DF{I);FOR  Jl-1  TO  I  DO 

BEGIN  U.HFU,J]i-M*H.HFlI'J]+MM*H.DF[I]*H.DF(J]  ; 
IF  I  O  J  THEN  U.HF[J,I]t-U.HF[I,J) 

END; 

END; 

END; 

RESl -U 
END; 
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APPENDIX  D 


Typical  output  of  the  program  SYSTEST. 

Values  computed  for  the  system  (6.4)  with  X[1)  -  X[2]  -  X[3) 

FUNCTION  VALUE: 

F[  1]  -  1 . 70000000000E+0 1 

FIRST  DERIVATIVES: 

DF I  1)/DX[  1)  -  6.40000000000E+01 

DF [  1]/DX[  2]  -  6. 4000000000 0E+0 1  ' 

DF [  11/DX(  3]  -  4 . OOOOOOOOOOOE+OO 

SECOND  DERIVATIVES: 

D2F  (  1)/DX(  1]DX(  1]  -  1 . 92 00000 0000E+02 

D2F [  1 ] /DX [  2 ] DX [  1]  >  0. OOOOOOOOOOOE+OO 

D2F t  1]/DX[  2 ] DX [  2]  -  1 .92000000000E+02 

D2F [  1]/DX[  3] DX [  1]  -  0. OOOOOOOOOOOE+OO 

D2F [  1]/DX[  3] DX [  2)  -  0. OOOOOOOOOOOE+OO 

D2F [  1 ] /DX [  3 ] DX [  3]  -  1 . 20000000000E+0 1 

FUNCTION  VALUE: 

F (  2]  -  0. OOOOOOOOOOOE+OO 

FIRST  DERIVATIVES: 

DF[  2]/DX[  1]  -  2. 00000000000E+00 

DF [  2]/DX[  2]  -  2. OOOOOOOOOOOE+OO 

DF [  2]/DX[  3]  -  2. OOOOOOOOOOOE+OO 

SECOND  DERIVATIVES: 

D2F [  2] /DX  [  1 ] DX [  1]  -  2. OOOOOOOOOOOE+OO 

D2F (  2]/DX(  2]DX[  1]  -  0. 00O00000000E+00 

D2F (  2]/DX[  2] DX [  2]  -  2.00000000000E+00 

D2F[  2]/DX[  3] DX [  1]  -  0. OOOOOOOOOOOE+OO 

D2F [  2]/DX[  3] DX [  2)  -  0. OOOOOOOOOOOE+OO 

D2F[  2 ] /DX [  3] DX [  3]  -  2. OOOOOOOOOOOE+OO 

FUNCTION  VALUE: 

Ft  3]  -  0. OOOOOOOOOOOE+OO 

FIRST  DERIVATIVES: 

DF [  3J/DX[  1]  -  3. OOOOOOOOOOOE+OO 

DF [  3] /DX [  2]  -  -1. OOOOOOOOOOOE+OO 
DF [  3] /DX [  3]  -  0. OOOOOOOOOOOE+OO 

SECOND  DERIVATIVES: 

D2F [  3] /DX [  1 ] DX [  1)  -  6. OOOOOOOOOOOE+OO 

D2F[  3]/DX[  2 ] DX [  1)  -  0. OOOOOOOOOOOE+OO 

D2F [  3]/DX[  2) DX [  2]  -  0. OOOOOOOOOOOE+OO 

D2F[  3]/DX(  3] DX [  1)  -  0. OOOOOOOOOOOE+OO 

D2F[  3]/DX[  3] DX [  2)  -  0. OOOOOOOOOOOE+OO 

D2F [  3]/DX(  3} DX [  3]  -  0. OOOOOOOOOOOE+OO 
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APPENDIX  E 


The  Paecal-SC  program  HALSYS  for  the  solution  of 
systems  of  equations  by  Halley's  method. 


Source  code  for  this  program  is  contained  in  the  file  HALSYS. PROG. 

PROGRAM  HALSYS  (INPUT,  OUTPUT)  I 

CONST  DIM  «  #i  (*  Replace  by  the  dimension  of  the  system  to  be  solved.  *) 

TYPE  DIMTYPE  “  1. .DIM; 

RVECTOR  -  ARRAY  [DIMTYPE]  OF  REAL ! 

RMATRIX  -  ARRAY  [DIMTYPE] OF  RVECTORi 

HESSIAN  -  RECORD  P ! REAL ) DF !  RVECTOR I HF I  RMATRIX  END I 
SYSTEM  -  ARRAY  [DIMTYPE]  OF  HESSIAN) 

VAR  It  DIMTYPElX, Ft  SYSTEM) Ct  CHAR,  A,B,Vt  RVECTOR) JAC,L, Ml  RMATRIX) 

NRSl BOOLEAN  I 


OPERATOR  +  (A,B:  RVECTOR)  RES:  RVECTOR) 

VAR  It  DIMTYPE; 

BEGIN  FOR  It-1  TO  DIM  DO  A[I]  :»  AtI]+B[I]| 

RES  1-  A 
END) 

OPERATOR  •  (At  REAL;  Bt  RVECTOR)  RES)  RVECTOR; 
VAR  It  DIMTYPE) 

BEGIN  FOR  It-1  TO  DIM  DO  B(I]  t»  A*B[U) 

RES  S 
END) 

OPERATOR  *  (At  RMATRIX)  Bt  RVECTOR)  RESt  RVECTORI 
VAR  It  DIMTYPE) 

BVARt  RVECTOR; 

BEGIN 
BVAR  1“  Bl 
FOR  I:*1  TO  DIM  DO 
BUI  I-  SCALP  (A  [I]  ,  BVAR,0  ) ; 

RES  !-  B 
END) 

OPERATOR  *  (A,Bt  RVECTOR)  RESt  RVECTORi 
VAR  It  DIMTYPE  I C I  RVECTOR) 

BEGIN  FOR  It-1  TO  DIM  DO  C[I]  t-A  [I]  *B[I)  » 
RISl-C 

END) 
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OPERATOR  /  (A,B:  RVECTOR)  RES:  RVECTOR) 

VAR  It  DIMTYPE |C:  RVECTOR; 

BEGIN  FOR  I:“  1  TO  DIM  DO 

IF  (A  [I]  "0 )  AND  (B(I]“0)  THEN  ClU:“0 
ELSE  ClI]»-AtI]/B[I); 

RES: “C 

END; 


FUNCTION  MRID:  RMATRIX;  <•  Returns  the  identity  matrix  *> 

VAR  I,J:  DIMTYPE; 

C:  RMATRIX; 

BEGIN 

FOR  I;“1  TO  DIM  DO 
FOR  J:«1  TO  DIM  DO 
IF  I«J  THEN  C(I,J]  1 
ELSE  C(I,Jj  0; 

MRID  :«  C 
END; 

FUNCTION  MRNULL:  RMATRIX;  (•  Returns  the  zero  matrix  •) 

VAR  I,J;  DIMTYPE; 

C:  RMATRIX; 

BEGIN 

FOR  Ii-I  TO  DIM  DO 
FOR  J;»1  TO  DIM  DO 
C[I,J)  0, 

MRNULL  C 
END; 

PROCEDURE  SOLVLN  <  D IM :  INTEGER ;  J  AC,  M:RMATRIX;V,S:  RVECTOR  ;NRS:  BOOLEAN); 

EXTERNAL  776; 

PROCEDURE  FCOMP(VAR  X.Fi  SYSTEM; DIM:  DIMTYPE); 

EXTERNAL  777; 

BEGIN  (*  Initialization  of  gradients  and  Hessians  of  independent  variables.  *) 

L: “MRID; FOR  I:-1  TO  DIM  DO 

BEGIN  Xll] .DF:“L[1J ;X [I] .HFi-MRNULL 

END; 

Ct“'R' iWHILE  C  “  'R'  DO 

BEGIN  (•  VALUE  INITIALIZATION  •) 

HRITELN CENTER  INDEPENDENT  VARIABLES'); 

FOR  I:“1  TO  DIM  DO  READ (X [1 1 . F ) ; 

Cl-'Q';  WHILE  C  -  'Q*  DO 
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BEGIN  (*  MAIN  PROGRAM  •) 


(•  Coapute  function  value*  and  darivativaai  print  value*  of  independent 
and  dapandant  varlablaa.  •) 

FCOMF(X,P,DIM)> 

FOR  Ii«1  TO  DIM  DO 

BEGIN  WRITEU((  'X(*(Ii2, ')  -  ’.XIU.F,*  F[',Xt2,M  -  ,,F(I].F)» 

END; 

WRITELN ('RESTART  (R)  OR  ITERATE  <Y/N)7’>) 

READ (C,C)i WHILE  C  *  '»'  SO 

BEGIN  <•  HALLEY  ITERATION  *) 

FOR  Il“1  TO  DIM  DO 

BEGIN  JACtI) t-F (I 1 .DF|  (•  JACOBIAN  MATRIX  •) 

V(U«— FCI].F|  (*  RIGHT  HAND  SIDE  •) 

END  I 

(*  Solve  for  A.  •) 

NR8t -FALSE) 80LVLM (DIM, JAC.M, V, A,NRS ) ) 

(*  Compute  B.  *) 

FOR  I«-1  TO  DIM  DO  L[U  l-Ftl)  .HF*A)V|-L«A) 

NRSi  -TRUElSOLVLN  (DIM,  JAC.M,  V,B,NRS  )  ) 

(•  Coapute  tha  Hallay  correction  and  update  independent  variables.  * ) 

V t  — A*A/ ( A*0 . 5*B) l FOR  I«-1  TO  JJIM  DO  X  [I  J  .F» -X  [I ]  .F+V [I ] ) 

Ci-'Q* 

END)  (•  HALLEY  ITERATION  •) 

END)  (*  MAIN  PROGRAM  •) 

END  (•  VALUE  INITIALIZATION  •) 

END. 
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APPENDIX  F 


Reaulta  of  the  Solution  of  the  System  (6.3)  by 
the  Program  HALSYS. 


INITIAL  VALUES 


xc 

1] 

- 

1.00000000000E+00 

PC 

U 

- 

1.70000000000E+01 

X[ 

2) 

m 

1.00000000000E+00 

Ft 

2] 

• 

O.OOOOOOOOOOOE+OO 

X[ 

3] 

m 

1.00000000000E+00 

Ft 

31 

“ 

0.00000000000E+00 

RESULTS 

OP  ITERATION  NUMBER  1 

X( 

1) 

a 

8.91 1 1870 1964E-0 1 

Ft 

1] 

a 

9 . 3752 1623 1 OOE-O 1 

X[ 

2] 

a 

7. 0 542934 1548E-01 

Ft 

2] 

» 

-9.44922445000E-03 

X[ 

3] 

" 

1.30339083879E+00 

Ft 

31 

■ 

2.20137281800E-03 

RESULTS 

OF  ITERATION  NUMBER  2 

X( 

1] 

m 

8. 77982528233E-0 1 

Ft 

11 

a 

1.036B5680000E-03 

X  ( 

2] 

- 

6.76786689302E-01 

Ft 

21 

- 

-9. 093240000 QOE-06 

x( 

3) 

m 

1.33082 58203 3E+00 

Ft 

3) 

“ 

9.05738500000E-06 

RESULTS 

OF  ITERATION  NUMBER  3 

x( 

1) 

a 

8.77965760275E-01 

Ft 

11 

a 

1.00000000000E-10 

X( 

2] 

- 

6. 767569705 19E-01 

Ft 

21 

a 

O.OOOOOOOOOOOE+OO 

xc 

31 

a 

1. 33085541 162E+00 

Ft 

31 

* 

O.OOOOOOOOOOOE+OO 

RESULTS 

OF  ITERATION  NUMBER  4 

XC 

1] 

m 

8. 779657602 74E-01 

rt 

11 

a 

O.OOOtlOOOOOOOE+OO 

X[ 

2] 

a 

6.767S6970516E-01 

Ft 

2] 

a 

O.OOOOOOOOOOOE+OO 

Xt 

31 

1. 33085541 162E+00 

Ft 

31 

a 

2.00000000000E-12 

RESULTS 

OF  ITERATION  NUMBER  5 

XI 

1] 

. 

8.77965  760274E-0 1 

Ft 

11 

a 

O.OOOOOOOOOOOE+OO 

x[ 

2) 

• 

6.76756970517E-01 

Ft 

21 

a 

O.OOOOOOOOOOOE+OO 

xc 

31 

a 

1. 33085541 162E+0O 

Ft 

31 

• 

1.00000000000E>12 

RESULTS 

OF  ITERATION  NUMBER  6 

X[ 

1) 

a 

8. 77965760 274E-01 

Ft 

11 

a 

O.OOOOOOOOOOOE+OO 

XI 

21 

• 

6. 767 569705 18E-01 

Ft 

2] 

- 

O.OOOOOOOOOOOE+OO 

XC 

31 

- 

1. 3308554 1162E+00 

Ft 

3! 

• 

O.OOOOOOOOOOOE+OO 
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